R Notebook: GLM and GLMM - Generalized Linear Models & Generalized Linear Mixed Models
Using Poisson, Negative Binomial and Linear Regression to Predict the Number of Visits and Lengths of Stay
Kyle P. Rasku RN BSN
07-15-2021
library(tidyverse)## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──## ✓ ggplot2 3.3.5 ✓ purrr 0.3.4
## ✓ tibble 3.1.2 ✓ dplyr 1.0.7
## ✓ tidyr 1.1.3 ✓ stringr 1.4.0
## ✓ readr 1.4.0 ✓ forcats 0.5.1## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()library(broom)
library(broom.mixed)
library(stats)
library(lme4)## Loading required package: Matrix##
## Attaching package: 'Matrix'## The following objects are masked from 'package:tidyr':
##
## expand, pack, unpacklibrary(lmerTest)##
## Attaching package: 'lmerTest'## The following object is masked from 'package:lme4':
##
## lmer## The following object is masked from 'package:stats':
##
## stepPUF_ELIX_IP2009 <- read.csv("PUF_ELIX_IP2009.csv")
PUF_ELIX_IP2009 <- within(PUF_ELIX_IP2009, {
STATE <- factor(STATE)
COUNTY <- factor(COUNTY)
FULL_FIPS_CODE <- factor(FULL_FIPS_CODE)
CBSA <- factor(CBSA)
})
PUF_ELIX_IP2009 <- PUF_ELIX_IP2009 %>%
mutate(VISITS = TOTAL_VISITS>0)BENE_AGE The beneficiary’s approximate age, calculated from the last day of the year 2009
B_SEX The beneficiary’s sex - 0 for males, 1 for females
B_DIED Whether or not this beneficiary died in 2009 (0 for No, 1 for Yes)
B_DIVERSE If the person is of color (1) or white (0)
HF_PROP_RANK The proportional rank of the person’s county against all the counties in their state on Health Factors (RWJF, 2010)
HO_PROP_RANK The proportional rank of the person’s county against all the counties in their state on Health Outcomes (RWJF, 2010)
TOTCHRONIC The max number of chronic conditions the beneficiary had in 2009
MEAN_ELIX_SCORE The mean Elixhauser score based on the person’s admissions in 2009 (Might be artificially 0 if a person with chronic conditions had no admissions in 2009.)
m1 <- glm(VISITS ~ BENE_AGE + B_SEX + B_DIVERSE + B_DIED + HF_PROP_RANK + HO_PROP_RANK + TOTCHRONIC + MEAN_ELIX_SCORE, data=PUF_ELIX_IP2009, family=binomial(link="logit"))
summary(m1)##
## Call:
## glm(formula = VISITS ~ BENE_AGE + B_SEX + B_DIVERSE + B_DIED +
## HF_PROP_RANK + HO_PROP_RANK + TOTCHRONIC + MEAN_ELIX_SCORE,
## family = binomial(link = "logit"), data = PUF_ELIX_IP2009)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.6333 -0.4555 -0.2246 -0.1804 3.0203
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.699990 0.067367 -54.923 < 2e-16 ***
## BENE_AGE -0.005143 0.000798 -6.445 1.16e-10 ***
## B_SEX 0.005323 0.020665 0.258 0.79674
## B_DIVERSE -0.051365 0.028395 -1.809 0.07046 .
## B_DIED 0.251596 0.078843 3.191 0.00142 **
## HF_PROP_RANK 0.057217 0.054299 1.054 0.29200
## HO_PROP_RANK -0.066118 0.057123 -1.157 0.24709
## TOTCHRONIC 0.413353 0.004435 93.200 < 2e-16 ***
## MEAN_ELIX_SCORE 0.475669 0.005182 91.794 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 93906 on 110329 degrees of freedom
## Residual deviance: 64097 on 110321 degrees of freedom
## (1852 observations deleted due to missingness)
## AIC: 64115
##
## Number of Fisher Scoring iterations: 6
In this model,
Age, Total # of Chronic Conditions and Mean Elixhauser Score are significant at the 99% CI against the outcome of having any inpatient “VISITS”.
As age goes up, the likelihood of any visit goes down. At some point, people die.
As the # of chronic conditions & mean Elixhauser scores increase the likelihood of any inpatient visits increase.
Whether or not the person Died in 2009 is significant at the 95% CI.
People often die in the hospital, so dying and inpatient visits are significantly related.
Not being white (B_DIVERSE) is significant at the 90% CI.
The estimate is negative, so there is a negative relationship between being not white and having visits.
Inverse relationships in correlation matrices show that there is a negative relationship between being a person of color and living longer, also system utilization / going to the doctor regularly.
Now that we’ve looked at the outcome of any “VISITS”, we can look at the outcome of TOTAL_VISITS>0
PUF_ELIX_IP2009_GT0 <- subset(PUF_ELIX_IP2009, TOTAL_VISITS>0)
The original data set had 113,305 observations, and this subset has 16,908.
First, let’s try OLS against the outcome of TOTAL_VISITS > 0. Then,
We will try a Poisson model with log link for the observations.
# Could be worse ... but standardized residual plots look pretty wonky!!
m0 <- lm(TOTAL_VISITS ~ BENE_AGE + B_SEX + B_DIVERSE + B_DIED + HF_PROP_RANK + HO_PROP_RANK + TOTCHRONIC + MEAN_ELIX_SCORE, data=PUF_ELIX_IP2009_GT0)
summary(m0)##
## Call:
## lm(formula = TOTAL_VISITS ~ BENE_AGE + B_SEX + B_DIVERSE + B_DIED +
## HF_PROP_RANK + HO_PROP_RANK + TOTCHRONIC + MEAN_ELIX_SCORE,
## data = PUF_ELIX_IP2009_GT0)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.0379 -0.3597 -0.1659 0.2019 2.8953
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.7239201 0.0287808 25.153 <2e-16 ***
## BENE_AGE -0.0002269 0.0003318 -0.684 0.494
## B_SEX 0.0069197 0.0085893 0.806 0.420
## B_DIVERSE 0.0189714 0.0118204 1.605 0.109
## B_DIED 0.3027413 0.0310598 9.747 <2e-16 ***
## HF_PROP_RANK 0.0006657 0.0224746 0.030 0.976
## HO_PROP_RANK -0.0140451 0.0235501 -0.596 0.551
## TOTCHRONIC 0.0870507 0.0020450 42.568 <2e-16 ***
## MEAN_ELIX_SCORE 0.0468164 0.0020604 22.722 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5467 on 16723 degrees of freedom
## (175 observations deleted due to missingness)
## Multiple R-squared: 0.1414, Adjusted R-squared: 0.141
## F-statistic: 344.2 on 8 and 16723 DF, p-value: < 2.2e-16plot(m0)
The OLS / Gaussian model is a poor fit.
m2 <- glm(TOTAL_VISITS ~ BENE_AGE + B_SEX + B_DIVERSE + B_DIED + HF_PROP_RANK + HO_PROP_RANK + TOTCHRONIC + MEAN_ELIX_SCORE, data=PUF_ELIX_IP2009_GT0, family=poisson(link = "log"))
summary(m2)##
## Call:
## glm(formula = TOTAL_VISITS ~ BENE_AGE + B_SEX + B_DIVERSE + B_DIED +
## HF_PROP_RANK + HO_PROP_RANK + TOTCHRONIC + MEAN_ELIX_SCORE,
## family = poisson(link = "log"), data = PUF_ELIX_IP2009_GT0)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.9085 -0.3100 -0.1467 0.1304 2.1249
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.2020799 0.0466156 -4.335 1.46e-05 ***
## BENE_AGE -0.0001485 0.0005293 -0.281 0.779
## B_SEX 0.0049213 0.0137592 0.358 0.721
## B_DIVERSE 0.0141917 0.0188282 0.754 0.451
## B_DIED 0.2043975 0.0447024 4.572 4.82e-06 ***
## HF_PROP_RANK 0.0002444 0.0359756 0.007 0.995
## HO_PROP_RANK -0.0099202 0.0376879 -0.263 0.792
## TOTCHRONIC 0.0681299 0.0033237 20.498 < 2e-16 ***
## MEAN_ELIX_SCORE 0.0357517 0.0032438 11.021 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 3669.4 on 16731 degrees of freedom
## Residual deviance: 3034.6 on 16723 degrees of freedom
## (175 observations deleted due to missingness)
## AIC: 39373
##
## Number of Fisher Scoring iterations: 4
Age is no longer significant in this model.
Death, Chronic Conditions and Elixhauser Score remain significant.
Let’s compare this Poisson model to a Negative Binomial model.
library(MASS)##
## Attaching package: 'MASS'## The following object is masked from 'package:dplyr':
##
## selectm3 <- glm.nb(TOTAL_VISITS ~ BENE_AGE + B_SEX + B_DIVERSE + B_DIED + HF_PROP_RANK + HO_PROP_RANK + TOTCHRONIC + MEAN_ELIX_SCORE, data=PUF_ELIX_IP2009_GT0)## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reachedsummary(m3)##
## Call:
## glm.nb(formula = TOTAL_VISITS ~ BENE_AGE + B_SEX + B_DIVERSE +
## B_DIED + HF_PROP_RANK + HO_PROP_RANK + TOTCHRONIC + MEAN_ELIX_SCORE,
## data = PUF_ELIX_IP2009_GT0, init.theta = 118176.1914, link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.9085 -0.3100 -0.1467 0.1304 2.1249
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.2020797 0.0466160 -4.335 1.46e-05 ***
## BENE_AGE -0.0001485 0.0005293 -0.281 0.779
## B_SEX 0.0049212 0.0137593 0.358 0.721
## B_DIVERSE 0.0141917 0.0188283 0.754 0.451
## B_DIED 0.2043973 0.0447029 4.572 4.82e-06 ***
## HF_PROP_RANK 0.0002444 0.0359759 0.007 0.995
## HO_PROP_RANK -0.0099201 0.0376882 -0.263 0.792
## TOTCHRONIC 0.0681299 0.0033237 20.498 < 2e-16 ***
## MEAN_ELIX_SCORE 0.0357517 0.0032438 11.021 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(118176.2) family taken to be 1)
##
## Null deviance: 3669.3 on 16731 degrees of freedom
## Residual deviance: 3034.5 on 16723 degrees of freedom
## (175 observations deleted due to missingness)
## AIC: 39375
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 118176
## Std. Err.: 170345
## Warning while fitting theta: iteration limit reached
##
## 2 x log-likelihood: -39354.98And finally, let’s fit the original data to a Zero-Inflated Negative Binomial model
library(pscl)## Classes and Methods for R developed in the
## Political Science Computational Laboratory
## Department of Political Science
## Stanford University
## Simon Jackman
## hurdle and zeroinfl functions by Achim Zeileism4 <- zeroinfl(TOTAL_VISITS ~ BENE_AGE + B_SEX + B_DIVERSE + B_DIED + HF_PROP_RANK + HO_PROP_RANK + TOTCHRONIC + MEAN_ELIX_SCORE,
data = PUF_ELIX_IP2009, dist = "negbin")
summary(m4)##
## Call:
## zeroinfl(formula = TOTAL_VISITS ~ BENE_AGE + B_SEX + B_DIVERSE + B_DIED +
## HF_PROP_RANK + HO_PROP_RANK + TOTCHRONIC + MEAN_ELIX_SCORE, data = PUF_ELIX_IP2009,
## dist = "negbin")
##
## Pearson residuals:
## Min 1Q Median 3Q Max
## -1.68015 -0.34172 -0.11074 -0.08544 30.89873
##
## Count model coefficients (negbin with log link):
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.2930409 0.0527404 -43.478 < 2e-16 ***
## BENE_AGE -0.0018493 0.0005757 -3.212 0.00132 **
## B_SEX 0.0246383 0.0150825 1.634 0.10235
## B_DIVERSE 0.0212150 0.0206331 1.028 0.30385
## B_DIED 0.3635194 0.0489567 7.425 1.13e-13 ***
## HF_PROP_RANK 0.0378638 0.0394239 0.960 0.33684
## HO_PROP_RANK -0.0648401 0.0414014 -1.566 0.11732
## TOTCHRONIC 0.2441284 0.0035428 68.908 < 2e-16 ***
## MEAN_ELIX_SCORE 0.1583228 0.0041000 38.615 < 2e-16 ***
## Log(theta) 14.0213433 18.6855791 0.750 0.45302
##
## Zero-inflation model coefficients (binomial with logit link):
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.591730 0.212105 7.504 6.17e-14 ***
## BENE_AGE 0.008951 0.002618 3.420 0.000627 ***
## B_SEX 0.102321 0.064578 1.584 0.113090
## B_DIVERSE 0.211687 0.087673 2.415 0.015756 *
## B_DIED 0.316981 0.241339 1.313 0.189038
## HF_PROP_RANK 0.047620 0.169204 0.281 0.778377
## HO_PROP_RANK -0.153982 0.178558 -0.862 0.388485
## TOTCHRONIC -0.223664 0.015084 -14.828 < 2e-16 ***
## MEAN_ELIX_SCORE -2.191743 0.080080 -27.369 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Theta = 1228547.7199
## Number of iterations in BFGS optimization: 34
## Log-likelihood: -4.173e+04 on 19 Df
As you can see, we could have jumped straight to the Zero-Inflated model, as it contains the hurdle model!
Let’s look at some of the plots.
# Diagnostic and Residuals plot for m2: Poisson regression (over dispersion is evident)
plot(m2)
# Diagnostic and Residuals plot for m3: Negative Binomial regression (still over dispersed)
plot(m3)
# Negative Binomial
ggplot(data = PUF_ELIX_IP2009_GT0, aes(x = TOTCHRONIC, y = TOTAL_VISITS)) +
geom_jitter(height=0.05, width = 0.1) +
stat_smooth(method = "glm.nb")## `geom_smooth()` using formula 'y ~ x'## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
# Poisson
ggplot(data = PUF_ELIX_IP2009_GT0, aes(x = TOTCHRONIC, y = TOTAL_VISITS)) +
geom_jitter(height=0.05, width = 0.1) +
stat_smooth(method = "glm", method.args = c(family="poisson"))## `geom_smooth()` using formula 'y ~ x'
They’re basically identical.
# Negative Binomial
ggplot(data = PUF_ELIX_IP2009_GT0, aes(x = MEAN_ELIX_SCORE, y = TOTAL_VISITS)) +
geom_jitter(height=0.05, width = 0.1) +
stat_smooth(method = "glm.nb")## `geom_smooth()` using formula 'y ~ x'## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
# Poisson
ggplot(data = PUF_ELIX_IP2009_GT0, aes(x = MEAN_ELIX_SCORE, y = TOTAL_VISITS)) +
geom_jitter(height=0.05, width = 0.1) +
stat_smooth(method = "glm", method.args = c(family="poisson"))## `geom_smooth()` using formula 'y ~ x'
These are very similar.
Based on the dispersion handling, the zero-inflated model was the most accurate thus far,
but let’s look at some different covariates, using the Poisson model.
What else could be related to the total # of inpatient visits?
# Can't use TOTAL_PHYS or TOTAL_CODES due to their near-perfect correlation with the outcome variable (see python corr matrix)
m2b <- glm(TOTAL_VISITS ~ BENE_AGE + B_DIVERSE + B_DIED + HF_PROP_RANK + HO_PROP_RANK + TOTCHRONIC + MEAN_ELIX_SCORE + CVRG_MOS + RX_CVRG_MOS + TOTAL_LOS, data=PUF_ELIX_IP2009_GT0, family=poisson(link = "log"))
summary(m2b)##
## Call:
## glm(formula = TOTAL_VISITS ~ BENE_AGE + B_DIVERSE + B_DIED +
## HF_PROP_RANK + HO_PROP_RANK + TOTCHRONIC + MEAN_ELIX_SCORE +
## CVRG_MOS + RX_CVRG_MOS + TOTAL_LOS, family = poisson(link = "log"),
## data = PUF_ELIX_IP2009_GT0)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.47013 -0.22946 -0.09558 0.07130 1.95714
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.654e-01 9.852e-02 -2.694 0.00706 **
## BENE_AGE -3.128e-05 5.262e-04 -0.059 0.95260
## B_DIVERSE 4.418e-03 1.884e-02 0.235 0.81455
## B_DIED 1.462e-01 4.476e-02 3.267 0.00109 **
## HF_PROP_RANK -1.167e-03 3.600e-02 -0.032 0.97414
## HO_PROP_RANK -3.371e-03 3.774e-02 -0.089 0.92883
## TOTCHRONIC 4.371e-02 3.450e-03 12.670 < 2e-16 ***
## MEAN_ELIX_SCORE 1.959e-02 3.317e-03 5.906 3.51e-09 ***
## CVRG_MOS 1.675e-03 7.327e-03 0.229 0.81915
## RX_CVRG_MOS 1.290e-03 1.720e-03 0.750 0.45322
## TOTAL_LOS 2.829e-02 8.923e-04 31.703 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 3669.4 on 16731 degrees of freedom
## Residual deviance: 2156.8 on 16721 degrees of freedom
## (175 observations deleted due to missingness)
## AIC: 38499
##
## Number of Fisher Scoring iterations: 4
Now, let’s see what impact State and County or CBSA would have.
m2c <- glm(TOTAL_VISITS ~ BENE_AGE + B_DIVERSE + B_DIED + HF_PROP_RANK + HO_PROP_RANK + TOTCHRONIC + MEAN_ELIX_SCORE + TOTAL_LOS + STATE, data=PUF_ELIX_IP2009_GT0, family=poisson(link = "log"))
summary(m2c)##
## Call:
## glm(formula = TOTAL_VISITS ~ BENE_AGE + B_DIVERSE + B_DIED +
## HF_PROP_RANK + HO_PROP_RANK + TOTCHRONIC + MEAN_ELIX_SCORE +
## TOTAL_LOS + STATE, family = poisson(link = "log"), data = PUF_ELIX_IP2009_GT0)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.45431 -0.22923 -0.09585 0.07105 1.95996
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.229e-01 6.572e-02 -3.391 0.000697 ***
## BENE_AGE -4.833e-05 5.282e-04 -0.091 0.927104
## B_DIVERSE 3.524e-03 1.929e-02 0.183 0.855078
## B_DIED 1.458e-01 4.489e-02 3.247 0.001165 **
## HF_PROP_RANK 3.099e-03 3.833e-02 0.081 0.935568
## HO_PROP_RANK -7.839e-03 4.055e-02 -0.193 0.846697
## TOTCHRONIC 4.386e-02 3.437e-03 12.759 < 2e-16 ***
## MEAN_ELIX_SCORE 1.957e-02 3.322e-03 5.892 3.83e-09 ***
## TOTAL_LOS 2.832e-02 8.951e-04 31.642 < 2e-16 ***
## STATEAlaska 3.373e-02 1.783e-01 0.189 0.849975
## STATEArizona -3.518e-02 7.256e-02 -0.485 0.627764
## STATEArkansas -3.369e-03 7.795e-02 -0.043 0.965530
## STATECalifornia -1.741e-03 5.368e-02 -0.032 0.974129
## STATEColorado -1.584e-02 7.701e-02 -0.206 0.837059
## STATEConnecticut -7.701e-04 8.125e-02 -0.009 0.992437
## STATEDelaware -8.329e-03 1.266e-01 -0.066 0.947555
## STATEDistrict of Columbia -4.154e-04 1.610e-01 -0.003 0.997941
## STATEFlorida -8.418e-03 5.432e-02 -0.155 0.876838
## STATEGeorgia 1.598e-03 6.204e-02 0.026 0.979450
## STATEHawaii -2.357e-03 1.261e-01 -0.019 0.985087
## STATEIdaho 6.570e-03 1.124e-01 0.058 0.953389
## STATEIllinois -9.086e-03 5.785e-02 -0.157 0.875207
## STATEIndiana -1.101e-02 6.505e-02 -0.169 0.865647
## STATEIowa -1.995e-02 7.905e-02 -0.252 0.800749
## STATEKansas -6.933e-02 8.411e-02 -0.824 0.409816
## STATEKentucky -9.691e-03 6.953e-02 -0.139 0.889152
## STATELouisiana -1.039e-02 7.232e-02 -0.144 0.885782
## STATEMaine -7.726e-03 9.895e-02 -0.078 0.937767
## STATEMaryland 7.252e-03 6.953e-02 0.104 0.916924
## STATEMassachusetts -2.879e-02 6.451e-02 -0.446 0.655398
## STATEMichigan -6.699e-03 5.911e-02 -0.113 0.909776
## STATEMinnesota 9.549e-03 7.240e-02 0.132 0.895067
## STATEMississippi -2.220e-02 7.367e-02 -0.301 0.763129
## STATEMissouri -1.074e-02 6.586e-02 -0.163 0.870468
## STATEMontana -3.387e-03 1.322e-01 -0.026 0.979564
## STATENebraska -7.422e-02 9.288e-02 -0.799 0.424217
## STATENevada -5.976e-02 9.925e-02 -0.602 0.547103
## STATENew Hampshire -8.496e-02 1.093e-01 -0.778 0.436822
## STATENew Jersey 6.841e-05 6.161e-02 0.001 0.999114
## STATENew Mexico -1.578e-02 1.064e-01 -0.148 0.882135
## STATENew York -1.869e-02 5.575e-02 -0.335 0.737502
## STATENorth Carolina -1.092e-02 6.091e-02 -0.179 0.857768
## STATENorth Dakota 8.357e-02 1.305e-01 0.640 0.521923
## STATEOhio -4.039e-03 5.844e-02 -0.069 0.944893
## STATEOklahoma 3.336e-02 7.448e-02 0.448 0.654260
## STATEOregon -2.961e-02 8.336e-02 -0.355 0.722473
## STATEPennsylvania 2.958e-03 5.747e-02 0.051 0.958948
## STATERhode Island -4.776e-02 1.255e-01 -0.381 0.703450
## STATESouth Carolina -2.662e-03 6.984e-02 -0.038 0.969596
## STATESouth Dakota -4.225e-02 1.297e-01 -0.326 0.744577
## STATETennessee -4.676e-02 6.523e-02 -0.717 0.473507
## STATETexas 1.154e-02 5.500e-02 0.210 0.833866
## STATEUtah -8.078e-03 1.036e-01 -0.078 0.937859
## STATEVermont -5.687e-02 1.552e-01 -0.366 0.714037
## STATEVirginia -3.215e-02 6.444e-02 -0.499 0.617847
## STATEWashington -4.231e-02 6.926e-02 -0.611 0.541304
## STATEWest Virginia 7.322e-02 8.493e-02 0.862 0.388607
## STATEWisconsin -5.365e-04 6.801e-02 -0.008 0.993706
## STATEWyoming -1.627e-02 1.568e-01 -0.104 0.917394
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 3669.4 on 16731 degrees of freedom
## Residual deviance: 2148.5 on 16673 degrees of freedom
## (175 observations deleted due to missingness)
## AIC: 38587
##
## Number of Fisher Scoring iterations: 4
No significant states. Let’s try CBSAs.
m2e <- glm(TOTAL_VISITS ~ BENE_AGE + B_DIVERSE + B_DIED + HF_PROP_RANK + HO_PROP_RANK + TOTCHRONIC + MEAN_ELIX_SCORE + RX_CVRG_MOS + TOTAL_LOS + CBSA, data=PUF_ELIX_IP2009_GT0, family=Gamma(link = "log"))
summary(m2e)##
## Call:
## glm(formula = TOTAL_VISITS ~ BENE_AGE + B_DIVERSE + B_DIED +
## HF_PROP_RANK + HO_PROP_RANK + TOTCHRONIC + MEAN_ELIX_SCORE +
## RX_CVRG_MOS + TOTAL_LOS + CBSA, family = Gamma(link = "log"),
## data = PUF_ELIX_IP2009_GT0)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.79865 -0.19922 -0.08349 0.06562 1.40956
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.200e-01 1.796e-02 -12.251 < 2e-16 ***
## BENE_AGE 1.427e-05 1.915e-04 0.075 0.94058
## B_DIVERSE 4.897e-03 6.934e-03 0.706 0.48004
## B_DIED 1.374e-01 1.793e-02 7.662 1.93e-14 ***
## HF_PROP_RANK 1.245e-02 1.840e-02 0.677 0.49865
## HO_PROP_RANK -1.783e-02 1.865e-02 -0.957 0.33882
## TOTCHRONIC 3.665e-02 1.225e-03 29.930 < 2e-16 ***
## MEAN_ELIX_SCORE 1.704e-02 1.207e-03 14.116 < 2e-16 ***
## RX_CVRG_MOS 8.316e-04 6.053e-04 1.374 0.16948
## TOTAL_LOS 3.161e-02 4.130e-04 76.528 < 2e-16 ***
## CBSA10180 2.063e-01 8.716e-02 2.367 0.01793 *
## CBSA10420 1.350e-02 4.864e-02 0.277 0.78141
## CBSA10500 -4.525e-02 9.031e-02 -0.501 0.61633
## CBSA10580 -3.526e-03 4.807e-02 -0.073 0.94153
## CBSA10740 5.629e-03 5.110e-02 0.110 0.91230
## CBSA10780 2.706e-02 9.439e-02 0.287 0.77438
## CBSA10900 -4.473e-02 4.071e-02 -1.099 0.27181
## CBSA11020 6.080e-02 9.922e-02 0.613 0.54002
## CBSA11100 -1.261e-01 8.681e-02 -1.452 0.14649
## CBSA11180 -1.334e-01 1.399e-01 -0.954 0.34026
## CBSA11260 3.009e-02 7.622e-02 0.395 0.69298
## CBSA11300 -8.264e-02 8.688e-02 -0.951 0.34147
## CBSA11340 -5.405e-02 8.690e-02 -0.622 0.53396
## CBSA11460 9.511e-02 7.404e-02 1.285 0.19897
## CBSA11500 -1.443e-02 9.900e-02 -0.146 0.88408
## CBSA11540 1.259e-02 1.398e-01 0.090 0.92826
## CBSA11700 -6.832e-02 5.497e-02 -1.243 0.21397
## CBSA12020 8.107e-02 1.399e-01 0.580 0.56222
## CBSA12060 5.941e-03 2.287e-02 0.260 0.79499
## CBSA12100 1.355e-02 9.042e-02 0.150 0.88087
## CBSA12220 4.768e-03 1.805e-01 0.026 0.97893
## CBSA12260 2.515e-03 5.164e-02 0.049 0.96116
## CBSA12420 1.120e-01 4.208e-02 2.661 0.00780 **
## CBSA12540 -5.175e-02 6.163e-02 -0.840 0.40110
## CBSA12580 3.042e-04 2.642e-02 0.012 0.99082
## CBSA12620 -1.091e-01 9.034e-02 -1.208 0.22720
## CBSA12700 3.685e-02 6.683e-02 0.551 0.58131
## CBSA12940 9.324e-03 5.406e-02 0.172 0.86307
## CBSA12980 1.344e-01 1.106e-01 1.215 0.22440
## CBSA13020 -2.486e-03 1.276e-01 -0.019 0.98446
## CBSA13140 -2.607e-02 6.156e-02 -0.423 0.67196
## CBSA13380 -1.036e-01 9.046e-02 -1.145 0.25228
## CBSA13460 7.974e-02 1.044e-01 0.764 0.44496
## CBSA13644 4.349e-02 4.572e-02 0.951 0.34146
## CBSA13740 -1.727e-02 9.895e-02 -0.175 0.86143
## CBSA13780 4.656e-02 7.193e-02 0.647 0.51748
## CBSA13820 1.906e-02 3.535e-02 0.539 0.58968
## CBSA13900 6.261e-03 9.901e-02 0.063 0.94958
## CBSA13980 -6.874e-02 8.087e-02 -0.850 0.39533
## CBSA14020 -1.059e-01 9.430e-02 -1.124 0.26123
## CBSA14060 1.149e-02 1.183e-01 0.097 0.92260
## CBSA14260 -4.319e-02 6.689e-02 -0.646 0.51847
## CBSA14484 -2.021e-02 3.071e-02 -0.658 0.51044
## CBSA14500 2.888e-02 9.050e-02 0.319 0.74961
## CBSA14540 7.729e-02 1.398e-01 0.553 0.58039
## CBSA14740 -1.168e-01 9.437e-02 -1.237 0.21602
## CBSA14860 -4.225e-02 5.551e-02 -0.761 0.44665
## CBSA15180 2.988e-02 6.651e-02 0.449 0.65327
## CBSA15260 1.557e-01 1.563e-01 0.996 0.31928
## CBSA15380 6.870e-02 4.076e-02 1.686 0.09191 .
## CBSA15500 -8.876e-02 1.182e-01 -0.751 0.45263
## CBSA15540 -2.147e-02 1.182e-01 -0.182 0.85587
## CBSA15764 -2.191e-02 3.329e-02 -0.658 0.51046
## CBSA15804 4.231e-04 3.828e-02 0.011 0.99118
## CBSA15940 3.923e-03 5.925e-02 0.066 0.94721
## CBSA15980 2.933e-03 4.748e-02 0.062 0.95074
## CBSA16020 7.334e-03 8.369e-02 0.088 0.93017
## CBSA16180 2.742e-01 1.562e-01 1.755 0.07921 .
## CBSA16220 5.879e-02 1.105e-01 0.532 0.59483
## CBSA16300 -6.216e-02 6.537e-02 -0.951 0.34169
## CBSA16580 -8.000e-02 1.276e-01 -0.627 0.53079
## CBSA16620 6.249e-04 6.561e-02 0.010 0.99240
## CBSA16700 1.656e-01 5.183e-02 3.194 0.00140 **
## CBSA16740 -3.385e-02 3.427e-02 -0.988 0.32329
## CBSA16820 -8.233e-02 7.835e-02 -1.051 0.29337
## CBSA16860 -6.139e-02 5.332e-02 -1.151 0.24958
## CBSA16940 3.225e-02 1.398e-01 0.231 0.81748
## CBSA16974 -1.001e-02 1.621e-02 -0.617 0.53703
## CBSA17020 -9.520e-03 1.182e-01 -0.081 0.93579
## CBSA17140 -1.956e-02 2.827e-02 -0.692 0.48895
## CBSA17300 1.421e-02 8.379e-02 0.170 0.86532
## CBSA17420 1.277e-01 8.094e-02 1.578 0.11466
## CBSA17460 -2.554e-02 2.778e-02 -0.920 0.35780
## CBSA17660 1.720e-01 1.276e-01 1.348 0.17776
## CBSA17780 2.605e-01 2.209e-01 1.179 0.23836
## CBSA17820 3.677e-03 5.823e-02 0.063 0.94965
## CBSA17860 3.178e-01 1.107e-01 2.871 0.00409 **
## CBSA17900 -2.428e-02 5.329e-02 -0.456 0.64862
## CBSA17980 -1.088e-01 8.681e-02 -1.253 0.21014
## CBSA18020 1.016e-01 1.106e-01 0.919 0.35834
## CBSA18140 2.945e-02 3.284e-02 0.897 0.36986
## CBSA18580 6.387e-02 5.924e-02 1.078 0.28099
## CBSA18700 3.088e-02 3.124e-01 0.099 0.92126
## CBSA18880 3.411e-02 1.044e-01 0.327 0.74386
## CBSA19060 -3.655e-02 1.277e-01 -0.286 0.77464
## CBSA19124 3.672e-02 2.629e-02 1.397 0.16254
## CBSA19140 -1.809e-01 1.107e-01 -1.634 0.10224
## CBSA19180 -2.675e-01 1.277e-01 -2.096 0.03614 *
## CBSA19260 1.383e-01 9.040e-02 1.530 0.12592
## CBSA19340 1.692e-01 7.008e-02 2.415 0.01575 *
## CBSA19380 -1.943e-03 4.323e-02 -0.045 0.96415
## CBSA19460 -8.210e-02 8.377e-02 -0.980 0.32710
## CBSA19500 -6.733e-02 1.277e-01 -0.527 0.59792
## CBSA19660 -2.127e-02 4.791e-02 -0.444 0.65703
## CBSA19740 -5.435e-03 3.116e-02 -0.174 0.86153
## CBSA19780 -1.536e-02 5.567e-02 -0.276 0.78266
## CBSA19804 -4.410e-02 2.801e-02 -1.574 0.11545
## CBSA20020 6.476e-02 8.688e-02 0.745 0.45608
## CBSA20100 -1.171e-01 9.912e-02 -1.181 0.23763
## CBSA20220 5.485e-02 1.182e-01 0.464 0.64273
## CBSA20260 -1.287e-01 8.375e-02 -1.537 0.12426
## CBSA20500 -1.839e-02 6.170e-02 -0.298 0.76567
## CBSA20740 -1.299e-01 1.182e-01 -1.099 0.27167
## CBSA20764 1.729e-02 2.808e-02 0.616 0.53808
## CBSA20940 2.337e-01 1.107e-01 2.112 0.03472 *
## CBSA21060 1.673e-01 1.107e-01 1.512 0.13061
## CBSA21140 3.301e-02 1.048e-01 0.315 0.75272
## CBSA21300 -1.290e-01 1.563e-01 -0.825 0.40928
## CBSA21340 -6.132e-03 5.375e-02 -0.114 0.90916
## CBSA21500 4.269e-02 7.637e-02 0.559 0.57614
## CBSA21660 -1.452e-01 8.079e-02 -1.797 0.07241 .
## CBSA21780 -5.174e-02 7.029e-02 -0.736 0.46162
## CBSA21820 -2.118e-01 1.804e-01 -1.174 0.24044
## CBSA22020 -1.231e-01 1.044e-01 -1.179 0.23830
## CBSA22140 -6.267e-02 1.805e-01 -0.347 0.72849
## CBSA22180 7.547e-02 7.184e-02 1.051 0.29348
## CBSA22220 1.637e-02 6.703e-02 0.244 0.80706
## CBSA22380 -8.428e-02 1.106e-01 -0.762 0.44606
## CBSA22420 -3.014e-03 6.415e-02 -0.047 0.96252
## CBSA22500 -7.560e-03 8.684e-02 -0.087 0.93062
## CBSA22520 -1.068e-01 9.896e-02 -1.079 0.28066
## CBSA22540 4.562e-02 1.276e-01 0.357 0.72077
## CBSA22660 -3.477e-02 6.414e-02 -0.542 0.58778
## CBSA22744 -1.128e-02 3.268e-02 -0.345 0.72999
## CBSA22900 4.958e-02 6.842e-02 0.725 0.46866
## CBSA23060 8.054e-02 6.151e-02 1.309 0.19040
## CBSA23104 4.179e-03 3.533e-02 0.118 0.90584
## CBSA23420 4.717e-02 5.655e-02 0.834 0.40418
## CBSA23460 -7.451e-02 1.042e-01 -0.715 0.47471
## CBSA23540 1.074e-01 8.685e-02 1.237 0.21622
## CBSA23580 2.630e-01 9.448e-02 2.783 0.00539 **
## CBSA23844 2.916e-02 4.746e-02 0.614 0.53893
## CBSA24020 7.196e-02 8.364e-02 0.860 0.38960
## CBSA24140 6.842e-02 1.276e-01 0.536 0.59185
## CBSA24220 -1.592e-01 1.277e-01 -1.247 0.21248
## CBSA24300 7.302e-02 9.032e-02 0.808 0.41883
## CBSA24340 -2.053e-02 4.374e-02 -0.469 0.63882
## CBSA24500 1.412e-01 1.106e-01 1.277 0.20168
## CBSA24540 -8.998e-02 1.563e-01 -0.576 0.56479
## CBSA24580 -6.194e-02 9.034e-02 -0.686 0.49296
## CBSA24660 -7.857e-02 4.983e-02 -1.577 0.11486
## CBSA24780 1.647e-01 1.563e-01 1.054 0.29194
## CBSA24860 -1.914e-02 4.714e-02 -0.406 0.68475
## CBSA25060 -1.813e-01 6.838e-02 -2.652 0.00802 **
## CBSA25180 1.432e-01 6.538e-02 2.190 0.02852 *
## CBSA25260 -2.269e-01 3.124e-01 -0.726 0.46767
## CBSA25420 9.453e-03 5.390e-02 0.175 0.86080
## CBSA25500 9.008e-02 9.440e-02 0.954 0.33999
## CBSA25540 3.946e-02 3.796e-02 1.040 0.29855
## CBSA25620 -4.831e-02 8.685e-02 -0.556 0.57804
## CBSA25860 5.716e-02 5.388e-02 1.061 0.28871
## CBSA25980 -3.831e-01 3.125e-01 -1.226 0.22017
## CBSA26100 -3.510e-02 6.706e-02 -0.523 0.60067
## CBSA26180 -6.240e-04 4.663e-02 -0.013 0.98932
## CBSA26300 -4.481e-02 8.392e-02 -0.534 0.59343
## CBSA26380 -6.305e-03 8.091e-02 -0.078 0.93789
## CBSA26420 2.744e-02 2.275e-02 1.206 0.22789
## CBSA26580 -3.036e-02 5.949e-02 -0.510 0.60985
## CBSA26620 4.109e-02 6.706e-02 0.613 0.54003
## CBSA26820 -7.029e-02 1.398e-01 -0.503 0.61515
## CBSA26900 2.607e-02 3.136e-02 0.831 0.40579
## CBSA26980 3.622e-01 1.805e-01 2.007 0.04481 *
## CBSA27060 -1.592e-01 1.563e-01 -1.018 0.30847
## CBSA27100 -9.886e-02 1.106e-01 -0.894 0.37135
## CBSA27140 1.808e-02 5.661e-02 0.319 0.74937
## CBSA27180 -6.696e-02 1.043e-01 -0.642 0.52081
## CBSA27260 -3.862e-02 3.944e-02 -0.979 0.32751
## CBSA27340 -1.816e-01 1.278e-01 -1.421 0.15531
## CBSA27500 9.930e-03 1.276e-01 0.078 0.93799
## CBSA27620 1.876e-02 1.106e-01 0.170 0.86528
## CBSA27740 4.750e-02 8.377e-02 0.567 0.57075
## CBSA27780 -1.093e-01 7.394e-02 -1.478 0.13948
## CBSA27860 -1.800e-02 2.209e-01 -0.081 0.93505
## CBSA27900 -8.769e-02 1.182e-01 -0.742 0.45823
## CBSA28020 -2.265e-02 6.548e-02 -0.346 0.72940
## CBSA28100 6.461e-02 1.182e-01 0.547 0.58461
## CBSA28140 1.565e-03 2.962e-02 0.053 0.95787
## CBSA28420 3.091e-02 8.366e-02 0.370 0.71174
## CBSA28660 7.286e-02 7.604e-02 0.958 0.33799
## CBSA28700 2.148e-02 5.935e-02 0.362 0.71746
## CBSA28740 2.066e-01 1.105e-01 1.869 0.06166 .
## CBSA28940 -3.699e-02 4.993e-02 -0.741 0.45887
## CBSA29020 -1.136e-01 1.105e-01 -1.028 0.30416
## CBSA29100 4.661e-02 9.055e-02 0.515 0.60673
## CBSA29140 7.764e-02 1.043e-01 0.744 0.45673
## CBSA29180 5.036e-02 8.380e-02 0.601 0.54785
## CBSA29340 3.252e-01 1.563e-01 2.081 0.03743 *
## CBSA29404 6.370e-02 4.757e-02 1.339 0.18055
## CBSA29420 8.849e-03 8.689e-02 0.102 0.91889
## CBSA29460 3.109e-02 4.596e-02 0.676 0.49875
## CBSA29540 -8.085e-02 5.748e-02 -1.407 0.15955
## CBSA29620 -6.053e-02 6.273e-02 -0.965 0.33462
## CBSA29700 -3.211e-02 9.150e-02 -0.351 0.72565
## CBSA29740 7.462e-03 9.921e-02 0.075 0.94004
## CBSA29820 -3.746e-02 3.888e-02 -0.964 0.33522
## CBSA29940 -2.719e-01 1.399e-01 -1.943 0.05202 .
## CBSA30020 9.360e-02 1.107e-01 0.846 0.39774
## CBSA30140 2.702e-02 1.277e-01 0.212 0.83236
## CBSA30300 -2.240e-01 1.804e-01 -1.242 0.21438
## CBSA30340 2.327e-01 1.398e-01 1.665 0.09587 .
## CBSA30460 -3.149e-02 5.949e-02 -0.529 0.59664
## CBSA30620 -1.250e-01 1.276e-01 -0.979 0.32738
## CBSA30700 -2.356e-02 7.026e-02 -0.335 0.73739
## CBSA30780 4.368e-02 5.480e-02 0.797 0.42546
## CBSA30860 -1.408e-01 2.209e-01 -0.637 0.52408
## CBSA30980 -2.462e-02 7.831e-02 -0.314 0.75325
## CBSA31020 -1.641e-01 1.398e-01 -1.174 0.24056
## CBSA31084 1.542e-02 1.912e-02 0.807 0.41993
## CBSA31140 -2.318e-02 3.833e-02 -0.605 0.54527
## CBSA31180 9.661e-02 1.043e-01 0.926 0.35446
## CBSA31340 -7.304e-02 7.829e-02 -0.933 0.35086
## CBSA31420 2.118e-02 1.106e-01 0.192 0.84812
## CBSA31460 1.320e-01 1.398e-01 0.944 0.34503
## CBSA31540 -6.989e-02 5.846e-02 -1.196 0.23190
## CBSA31700 -8.889e-02 7.005e-02 -1.269 0.20447
## CBSA31740 -2.770e-01 1.804e-01 -1.536 0.12466
## CBSA31860 -4.094e-02 1.277e-01 -0.321 0.74845
## CBSA31900 -2.156e-01 1.105e-01 -1.951 0.05110 .
## CBSA32580 -7.456e-02 5.577e-02 -1.337 0.18124
## CBSA32780 2.421e-01 1.182e-01 2.049 0.04045 *
## CBSA32820 1.774e-02 3.939e-02 0.450 0.65246
## CBSA32900 4.059e-02 1.563e-01 0.260 0.79507
## CBSA33124 -3.579e-02 2.948e-02 -1.214 0.22474
## CBSA33140 -9.469e-02 1.398e-01 -0.677 0.49816
## CBSA33260 -2.298e-01 1.398e-01 -1.643 0.10042
## CBSA33340 1.244e-02 3.229e-02 0.385 0.70017
## CBSA33460 -7.118e-03 2.790e-02 -0.255 0.79859
## CBSA33540 -3.547e-02 1.183e-01 -0.300 0.76424
## CBSA33660 -3.726e-02 6.554e-02 -0.569 0.56970
## CBSA33700 9.812e-03 6.537e-02 0.150 0.88069
## CBSA33740 -1.647e-01 9.433e-02 -1.746 0.08077 .
## CBSA33780 3.138e-01 1.106e-01 2.838 0.00455 **
## CBSA33860 -6.062e-02 6.155e-02 -0.985 0.32470
## CBSA34060 -9.938e-02 1.278e-01 -0.778 0.43673
## CBSA34100 1.771e-02 9.042e-02 0.196 0.84470
## CBSA34580 1.653e-01 9.891e-02 1.671 0.09473 .
## CBSA34620 -8.670e-02 1.043e-01 -0.831 0.40576
## CBSA34740 -2.334e-02 8.369e-02 -0.279 0.78029
## CBSA34820 -7.887e-02 7.597e-02 -1.038 0.29921
## CBSA34900 -1.272e-02 3.123e-01 -0.041 0.96752
## CBSA34940 8.095e-02 5.415e-02 1.495 0.13496
## CBSA34980 -1.845e-02 3.331e-02 -0.554 0.57963
## CBSA35004 -3.200e-02 2.602e-02 -1.230 0.21886
## CBSA35084 -6.104e-03 2.809e-02 -0.217 0.82800
## CBSA35300 2.518e-03 4.750e-02 0.053 0.95773
## CBSA35380 -3.278e-02 4.079e-02 -0.804 0.42160
## CBSA35644 -1.822e-02 1.440e-02 -1.265 0.20591
## CBSA35660 1.230e-01 9.431e-02 1.304 0.19227
## CBSA35840 -4.103e-02 3.678e-02 -1.116 0.26452
## CBSA35980 -1.940e-03 7.824e-02 -0.025 0.98022
## CBSA36084 -3.908e-02 3.370e-02 -1.160 0.24624
## CBSA36100 -2.020e-02 6.272e-02 -0.322 0.74734
## CBSA36140 -7.288e-02 9.436e-02 -0.772 0.43990
## CBSA36220 -4.370e-02 1.182e-01 -0.370 0.71161
## CBSA36260 4.754e-02 7.384e-02 0.644 0.51973
## CBSA36420 5.650e-02 3.809e-02 1.483 0.13800
## CBSA36500 1.404e-01 9.440e-02 1.488 0.13690
## CBSA36540 -4.502e-02 4.687e-02 -0.961 0.33680
## CBSA36740 -2.664e-02 2.904e-02 -0.917 0.35901
## CBSA36780 -1.353e-01 9.454e-02 -1.431 0.15231
## CBSA36980 -8.655e-03 9.440e-02 -0.092 0.92695
## CBSA37100 -1.050e-02 5.104e-02 -0.206 0.83703
## CBSA37340 5.896e-02 4.262e-02 1.383 0.16659
## CBSA37380 1.774e-01 1.182e-01 1.501 0.13348
## CBSA37460 -1.550e-01 9.435e-02 -1.642 0.10052
## CBSA37620 1.127e-02 9.436e-02 0.119 0.90489
## CBSA37700 3.985e-02 9.437e-02 0.422 0.67285
## CBSA37764 -4.762e-02 4.455e-02 -1.069 0.28519
## CBSA37860 -1.434e-02 4.965e-02 -0.289 0.77279
## CBSA37900 9.432e-02 5.567e-02 1.694 0.09023 .
## CBSA37964 2.685e-02 2.278e-02 1.178 0.23869
## CBSA38060 -2.969e-02 2.604e-02 -1.140 0.25431
## CBSA38220 2.313e-01 1.563e-01 1.480 0.13888
## CBSA38300 2.492e-02 2.521e-02 0.989 0.32287
## CBSA38340 2.879e-03 9.898e-02 0.029 0.97680
## CBSA38540 -7.889e-02 1.277e-01 -0.618 0.53673
## CBSA38860 1.769e-02 5.329e-02 0.332 0.73995
## CBSA38900 -3.140e-02 3.183e-02 -0.986 0.32396
## CBSA38940 -6.914e-02 5.557e-02 -1.244 0.21340
## CBSA39100 8.094e-02 5.475e-02 1.478 0.13930
## CBSA39140 5.048e-03 8.702e-02 0.058 0.95374
## CBSA39300 2.622e-03 3.165e-02 0.083 0.93398
## CBSA39340 -1.796e-01 9.442e-02 -1.903 0.05712 .
## CBSA39380 -1.146e-01 1.106e-01 -1.037 0.29998
## CBSA39460 -1.190e-01 8.369e-02 -1.422 0.15517
## CBSA39540 -9.359e-02 8.088e-02 -1.157 0.24723
## CBSA39580 4.262e-02 5.482e-02 0.777 0.43694
## CBSA39660 -6.808e-03 9.431e-02 -0.072 0.94246
## CBSA39740 -2.828e-03 7.013e-02 -0.040 0.96784
## CBSA39820 3.652e-02 1.105e-01 0.330 0.74108
## CBSA39900 -9.573e-02 7.015e-02 -1.365 0.17242
## CBSA40060 -6.924e-02 4.002e-02 -1.730 0.08364 .
## CBSA40140 1.062e-02 2.831e-02 0.375 0.70757
## CBSA40220 -4.389e-02 7.003e-02 -0.627 0.53080
## CBSA40340 -1.509e-01 1.277e-01 -1.182 0.23722
## CBSA40380 3.795e-02 3.755e-02 1.011 0.31223
## CBSA40420 -5.278e-03 6.401e-02 -0.082 0.93429
## CBSA40484 -5.887e-02 6.036e-02 -0.975 0.32943
## CBSA40580 -9.393e-02 1.562e-01 -0.601 0.54772
## CBSA40660 -3.809e-02 9.896e-02 -0.385 0.70032
## CBSA40900 4.263e-02 3.393e-02 1.257 0.20889
## CBSA40980 -6.363e-02 9.440e-02 -0.674 0.50026
## CBSA41060 -3.540e-03 1.043e-01 -0.034 0.97292
## CBSA41100 7.611e-03 1.398e-01 0.054 0.95658
## CBSA41140 3.651e-02 1.043e-01 0.350 0.72635
## CBSA41180 -2.804e-02 2.482e-02 -1.130 0.25859
## CBSA41420 1.483e-01 9.042e-02 1.641 0.10092
## CBSA41500 1.055e-01 8.690e-02 1.215 0.22454
## CBSA41540 -1.738e-01 1.106e-01 -1.571 0.11610
## CBSA41620 6.831e-03 4.291e-02 0.159 0.87351
## CBSA41660 -1.954e-02 9.038e-02 -0.216 0.82883
## CBSA41700 -4.552e-03 3.325e-02 -0.137 0.89112
## CBSA41740 -4.266e-02 2.977e-02 -1.433 0.15186
## CBSA41780 6.907e-02 1.805e-01 0.383 0.70199
## CBSA41884 -4.202e-02 3.204e-02 -1.312 0.18966
## CBSA41940 7.531e-02 4.247e-02 1.774 0.07616 .
## CBSA42020 -8.674e-02 8.376e-02 -1.036 0.30041
## CBSA42044 1.263e-02 2.939e-02 0.430 0.66737
## CBSA42060 -2.051e-02 5.931e-02 -0.346 0.72948
## CBSA42100 -6.000e-02 9.904e-02 -0.606 0.54465
## CBSA42140 -1.877e-01 1.277e-01 -1.469 0.14175
## CBSA42220 3.460e-02 6.414e-02 0.539 0.58962
## CBSA42340 6.356e-02 7.197e-02 0.883 0.37716
## CBSA42540 -9.764e-02 4.624e-02 -2.112 0.03471 *
## CBSA42644 -3.523e-02 3.018e-02 -1.167 0.24311
## CBSA42680 -3.958e-02 9.045e-02 -0.438 0.66170
## CBSA43100 -1.072e-01 1.182e-01 -0.907 0.36460
## CBSA43300 -1.849e-01 1.563e-01 -1.183 0.23699
## CBSA43340 1.269e-02 6.154e-02 0.206 0.83669
## CBSA43580 1.943e-01 1.398e-01 1.390 0.16466
## CBSA43620 -6.896e-02 8.097e-02 -0.852 0.39439
## CBSA43780 8.605e-03 7.824e-02 0.110 0.91243
## CBSA43900 -6.890e-02 7.388e-02 -0.933 0.35100
## CBSA44060 3.548e-02 6.402e-02 0.554 0.57947
## CBSA44100 -4.968e-02 7.832e-02 -0.634 0.52589
## CBSA44140 -1.215e-02 4.182e-02 -0.291 0.77134
## CBSA44180 8.930e-03 6.169e-02 0.145 0.88491
## CBSA44220 2.007e-01 1.562e-01 1.285 0.19890
## CBSA44300 2.067e-02 1.044e-01 0.198 0.84306
## CBSA44600 2.170e-02 9.035e-02 0.240 0.81018
## CBSA44700 -5.007e-02 6.851e-02 -0.731 0.46495
## CBSA44940 -6.766e-02 1.043e-01 -0.649 0.51656
## CBSA45060 -1.889e-02 5.239e-02 -0.361 0.71841
## CBSA45104 -3.663e-02 5.093e-02 -0.719 0.47207
## CBSA45220 -5.072e-02 7.830e-02 -0.648 0.51716
## CBSA45300 4.195e-03 2.457e-02 0.171 0.86442
## CBSA45460 -1.536e-01 8.682e-02 -1.769 0.07694 .
## CBSA45500 7.446e-02 1.042e-01 0.714 0.47495
## CBSA45780 1.436e-02 5.029e-02 0.286 0.77521
## CBSA45820 1.336e-01 1.398e-01 0.956 0.33929
## CBSA45940 6.435e-02 7.828e-02 0.822 0.41108
## CBSA46060 -7.105e-02 4.226e-02 -1.681 0.09270 .
## CBSA46140 2.594e-03 4.692e-02 0.055 0.95590
## CBSA46220 -7.421e-02 7.391e-02 -1.004 0.31539
## CBSA46340 -3.413e-02 7.598e-02 -0.449 0.65334
## CBSA46540 -1.121e-02 8.702e-02 -0.129 0.89754
## CBSA46660 -5.514e-02 1.182e-01 -0.467 0.64082
## CBSA46700 -9.762e-02 7.825e-02 -1.248 0.21220
## CBSA47020 -2.068e-01 1.277e-01 -1.619 0.10552
## CBSA47220 -1.161e-01 1.183e-01 -0.982 0.32611
## CBSA47260 8.331e-03 3.370e-02 0.247 0.80477
## CBSA47300 -3.046e-02 9.043e-02 -0.337 0.73621
## CBSA47380 9.507e-02 9.043e-02 1.051 0.29313
## CBSA47580 -1.627e-01 1.183e-01 -1.376 0.16889
## CBSA47644 7.739e-03 2.665e-02 0.290 0.77150
## CBSA47894 -2.308e-02 2.555e-02 -0.903 0.36630
## CBSA47940 3.127e-02 1.043e-01 0.300 0.76431
## CBSA48140 -1.786e-01 1.398e-01 -1.277 0.20154
## CBSA48300 1.892e-02 1.804e-01 0.105 0.91650
## CBSA48424 3.792e-02 3.422e-02 1.108 0.26774
## CBSA48540 -1.080e-02 9.031e-02 -0.120 0.90479
## CBSA48620 -1.040e-01 5.478e-02 -1.898 0.05769 .
## CBSA48660 3.701e-02 7.878e-02 0.470 0.63849
## CBSA48700 -2.247e-02 1.182e-01 -0.190 0.84921
## CBSA48864 8.832e-03 5.233e-02 0.169 0.86597
## CBSA48900 6.705e-02 5.739e-02 1.168 0.24265
## CBSA49020 -1.154e-01 1.398e-01 -0.826 0.40909
## CBSA49180 -7.099e-02 6.843e-02 -1.037 0.29957
## CBSA49340 -7.242e-02 5.389e-02 -1.344 0.17900
## CBSA49420 -8.664e-02 9.894e-02 -0.876 0.38122
## CBSA49620 8.453e-02 7.389e-02 1.144 0.25266
## CBSA49660 -3.370e-02 5.313e-02 -0.634 0.52594
## CBSA49700 -1.931e-01 1.563e-01 -1.236 0.21658
## CBSA49740 2.856e-02 1.284e-01 0.222 0.82395
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Gamma family taken to be 0.09748577)
##
## Null deviance: 2446.5 on 16731 degrees of freedom
## Residual deviance: 1342.3 on 16338 degrees of freedom
## (175 observations deleted due to missingness)
## AIC: 12691
##
## Number of Fisher Scoring iterations: 5
CBSA’s significant for the outcome of TOTAL_VISITS>0 include:
| 99% CI | 12420 | Austin-Round Rock, TX | POSITIVE |
| 16700 | Charleston-North Charleston, SC | POSITIVE | |
| 17860 | Columbia, MO | POSITIVE | |
| 23580 | Gainesville, GA | POSITIVE | |
| 25060 | Gulfport-Biloxi-Pascagoula, MS | NEGATIVE | |
| 95% CI | 10180 | Abilene, TX | POSITIVE |
| 19180 | Danville, IL | NEGATIVE | |
| 19340 | Davenport-Moline-Rock Island, IA-IL | POSITIVE | |
| 20940 | El Centro, CA | POSITIVE | |
| 25180 | Hagerstown-Martinsburg, MD-WV | POSITIVE | |
| 26980 | Iowa City, IA | POSITIVE | |
| 29340 | Lake Charles, LA | POSITIVE | |
| 90% CI | 15380 | Buffalo-Cheektowaga-Niagara Falls, NY | POSITIVE |
| 16180 | Carson City, NV | POSITIVE | |
| 21660 | Eugene, OR | NEGATIVE | |
| 28740 | Kingston, NY | POSITIVE | |
Now, let’s plot the outcome of TOTAL_VISITS against the MEAN_ELIX_SCORE, creating a separate intercept for each CBSA.
ggplot(data=PUF_ELIX_IP2009_GT0, aes(x=MEAN_ELIX_SCORE, y=TOTAL_VISITS,
group=CBSA)) +
geom_line() +
geom_smooth(method = "glm", method.args=c(family="poisson"), se=FALSE)## `geom_smooth()` using formula 'y ~ x'
#
# Poisson against the outcome of TOTAL_LOS (Total Length of Stay) for each CBSA
ggplot(data=PUF_ELIX_IP2009_GT0, aes(x=MEAN_ELIX_SCORE, y=TOTAL_LOS,
group=CBSA)) +
geom_line() +
geom_smooth(method = "glm", method.args=c(family="poisson"), se=FALSE)## `geom_smooth()` using formula 'y ~ x'
Now, let’s try fitting some models with the random effect of CBSA (noting that there are several beneficiaries whose counties are too small to belong to a CBSA)
zero_cbsas <- PUF_ELIX_IP2009 %>% subset(CBSA=="0")
zero_fips <- PUF_ELIX_IP2009 %>% subset(FULL_FIPS_CODE=="0")
25,943 beneficiaries without CBSAs.
1,755 beneficiaries without FIPS.
PUF_ELIX_IP2009_GT0 <- PUF_ELIX_IP2009_GT0 %>%
mutate(MC_SCALED = scale(MEAN_CHRONIC))
PUF_ELIX_IP2009_GT0 <- PUF_ELIX_IP2009_GT0 %>%
mutate(TC_SCALED = scale(TOTAL_CODES))
PUF_ELIX_IP2009_GT0 <- PUF_ELIX_IP2009_GT0 %>%
mutate(LOS_SCALED = scale(TOTAL_LOS))
PUF_ELIX_IP2009_GT0 <- PUF_ELIX_IP2009_GT0 %>%
mutate(PHYS_SCALED = scale(TOTAL_PHYS))
PUF_ELIX_IP2009_GT0 <- PUF_ELIX_IP2009_GT0 %>%
mutate(ELIX_SCALED = scale(MEAN_ELIX_SCORE))
CBSA_GT0 <- subset(PUF_ELIX_IP2009_GT0, CBSA!="0")
Modeling 13,080 Beneficiaries with 1 or more inpatient visits in 2009 who live in a CBSA:
# Linear Random Effects Model ~ RE Slope: Elixhauser Score, RE Intercept: CBSA
lmm = lmer(TOTAL_VISITS ~ MC_SCALED + LOS_SCALED + ELIX_SCALED + B_DIED + (ELIX_SCALED|CBSA), data = CBSA_GT0) ## boundary (singular) fit: see ?isSingularsummary(lmm)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: TOTAL_VISITS ~ MC_SCALED + LOS_SCALED + ELIX_SCALED + B_DIED +
## (ELIX_SCALED | CBSA)
## Data: CBSA_GT0
##
## REML criterion at convergence: 16922
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.6815 -0.5634 -0.1767 0.2737 5.8067
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## CBSA (Intercept) 0.0002498 0.01580
## ELIX_SCALED 0.0002179 0.01476 1.00
## Residual 0.2124769 0.46095
## Number of obs: 13080, groups: CBSA, 384
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 1.297e+00 4.292e-03 1.824e+02 302.116 < 2e-16 ***
## MC_SCALED 1.104e-01 4.192e-03 1.307e+04 26.331 < 2e-16 ***
## LOS_SCALED 3.032e-01 4.222e-03 1.307e+04 71.807 < 2e-16 ***
## ELIX_SCALED 5.054e-02 4.330e-03 2.062e+02 11.672 < 2e-16 ***
## B_DIED 2.188e-01 2.921e-02 1.307e+04 7.491 7.27e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) MC_SCA LOS_SC ELIX_S
## MC_SCALED -0.001
## LOS_SCALED 0.004 -0.232
## ELIX_SCALED 0.095 -0.087 -0.156
## B_DIED -0.132 0.002 -0.045 0.001
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingularplot(lmm)
# GLM Random Effects Model, Gamma ~ RE Intercept: CBSA
glmm = glmer(TOTAL_VISITS ~ MC_SCALED + LOS_SCALED + ELIX_SCALED + B_DIED + (1|CBSA), data = CBSA_GT0, family = Gamma(link = "log"))
summary(glmm)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: Gamma ( log )
## Formula: TOTAL_VISITS ~ MC_SCALED + LOS_SCALED + ELIX_SCALED + B_DIED +
## (1 | CBSA)
## Data: CBSA_GT0
##
## AIC BIC logLik deviance df.resid
## 9713.2 9765.6 -4849.6 9699.2 13073
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.8196 -0.5961 -0.2683 0.1865 6.5000
##
## Random effects:
## Groups Name Variance Std.Dev.
## CBSA (Intercept) 0.001422 0.03771
## Residual 0.098513 0.31387
## Number of obs: 13080, groups: CBSA, 384
##
## Fixed effects:
## Estimate Std. Error t value Pr(>|z|)
## (Intercept) 0.227911 0.004398 51.825 < 2e-16 ***
## MC_SCALED 0.077905 0.002604 29.915 < 2e-16 ***
## LOS_SCALED 0.193077 0.002766 69.793 < 2e-16 ***
## ELIX_SCALED 0.037110 0.002605 14.243 < 2e-16 ***
## B_DIED 0.132857 0.018120 7.332 2.27e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) MC_SCA LOS_SC ELIX_S
## MC_SCALED 0.002
## LOS_SCALED 0.003 -0.263
## ELIX_SCALED 0.004 -0.104 -0.178
## B_DIED -0.079 -0.003 -0.046 -0.002plot(glmm)
anova(lmm, glmm)## refitting model(s) with ML (instead of REML)In this case, the Gamma mixed model has the best fit, with AIC 9713.2 and -log-likelihood -4849.6
Both models have the same number of parameters.
CBSA_GT0 <- CBSA_GT0 %>%
mutate(SCALED_COSTS = scale(ALLCOSTS))
# Random Effects Model ~ RE Slope: TOTAL_PHYS (Number of Physicians), RE Intercept: CBSA
lmm2 = lmer(TOTAL_LOS ~ MEAN_CHRONIC + TOTAL_CODES + TOTAL_PHYS + MEAN_ELIX_SCORE + SCALED_COSTS + (TOTAL_PHYS|CBSA), data = CBSA_GT0) ## boundary (singular) fit: see ?isSingularsummary(lmm2)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula:
## TOTAL_LOS ~ MEAN_CHRONIC + TOTAL_CODES + TOTAL_PHYS + MEAN_ELIX_SCORE +
## SCALED_COSTS + (TOTAL_PHYS | CBSA)
## Data: CBSA_GT0
##
## REML criterion at convergence: 77595.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.4172 -0.5529 -0.1889 0.3470 9.0341
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## CBSA (Intercept) 4.1445 2.0358
## TOTAL_PHYS 0.9028 0.9502 -1.00
## Residual 21.4792 4.6346
## Number of obs: 13080, groups: CBSA, 384
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 9.136e-01 2.170e-01 6.463e+02 4.211 2.91e-05 ***
## MEAN_CHRONIC 1.172e-01 2.097e-02 1.287e+04 5.590 2.31e-08 ***
## TOTAL_CODES 3.303e-01 1.472e-02 1.305e+04 22.447 < 2e-16 ***
## TOTAL_PHYS 2.405e-01 9.904e-02 8.459e+02 2.428 0.0154 *
## MEAN_ELIX_SCORE 9.563e-02 2.058e-02 1.287e+04 4.646 3.42e-06 ***
## SCALED_COSTS 1.717e+00 7.058e-02 1.284e+04 24.330 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) MEAN_C TOTAL_C TOTAL_P MEAN_E
## MEAN_CHRONI -0.310
## TOTAL_CODES -0.064 -0.160
## TOTAL_PHYS -0.659 -0.034 -0.535
## MEAN_ELIX_S -0.181 -0.041 -0.180 0.048
## SCALED_COST 0.516 0.057 -0.328 -0.233 0.023
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingularplot(lmm2)
CBSA plus the slope of number of physicians explains the majority of the variance.
The remaining variance is small, but is well-explained by the covariates.
The residuals are fairly equally distributed around the zero line.
Let’s scale variables, and compare it with a glmer.
# GLM Random Effects Model ~ RE Slope: PHYS_SCALED (Number of Physicians), RE Intercept: CBSA
glmm2 = glmer.nb(TOTAL_LOS ~ MC_SCALED + TC_SCALED + PHYS_SCALED + ELIX_SCALED + SCALED_COSTS + (PHYS_SCALED|CBSA), data = CBSA_GT0) ## boundary (singular) fit: see ?isSingularsummary(glmm2)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: Negative Binomial(3.4561) ( log )
## Formula: TOTAL_LOS ~ MC_SCALED + TC_SCALED + PHYS_SCALED + ELIX_SCALED +
## SCALED_COSTS + (PHYS_SCALED | CBSA)
## Data: CBSA_GT0
##
## AIC BIC logLik deviance df.resid
## 69576.5 69651.3 -34778.3 69556.5 13070
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.6393 -0.7309 -0.2586 0.4263 8.5202
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## CBSA (Intercept) 7.041e-04 2.654e-02
## PHYS_SCALED 1.698e-10 1.303e-05 -0.87
## Number of obs: 13080, groups: CBSA, 384
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.743218 0.006475 269.202 < 2e-16 ***
## MC_SCALED 0.049163 0.006550 7.505 6.12e-14 ***
## TC_SCALED 0.332325 0.013244 25.093 < 2e-16 ***
## PHYS_SCALED 0.003009 0.012166 0.247 0.805
## ELIX_SCALED 0.040299 0.006256 6.442 1.18e-10 ***
## SCALED_COSTS 0.121043 0.009575 12.641 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) MC_SCA TC_SCA PHYS_S ELIX_S
## MC_SCALED -0.019
## TC_SCALED -0.052 -0.156
## PHYS_SCALED 0.006 -0.054 -0.677
## ELIX_SCALED -0.013 -0.025 -0.182 0.059
## SCALED_COST -0.005 0.056 -0.355 -0.262 0.024
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingularplot(glmm2)
ANOVA Comparison:
anova(lmm2, glmm2)## refitting model(s) with ML (instead of REML)
The glmer NB model has the superior fit, but this model doesn’t fit nearly as well as the visits model.
# LA COUNTY 06 037
county_subset <- CBSA_GT0 %>%
mutate(predictions = predict(glmm)) %>%
filter(FULL_FIPS_CODE=="6037") Exploring LA County’s 320 Beneficiaries, using the TOTAL_VISITS Gamma model
ggplot(county_subset, aes(x = predictions, y = TOTAL_VISITS)) +
geom_point() +
xlab("Predicted TOTAL VISITS (glmer Gamma)") +
ylab("Total Visits") +
theme_classic()
ggplot(county_subset, aes(x =predictions, y = (TOTAL_VISITS-predictions))) +
geom_point() +
xlab("Predicted Total Visits (glmer Gamma)") +
ylab("Residuals") +
theme_classic()
The best-fitting model still under-predicts in general
ggplot(county_subset, aes(x = MC_SCALED, y = TOTAL_VISITS)) +
geom_point() +
stat_smooth(method = "lm", aes(color="linear"), se=FALSE) +
stat_smooth(method = "glm", aes(color="gamma"), method.args=list(family = Gamma(link = "log")), se=FALSE) +
stat_smooth(method = "glm.nb", aes(color="negative binomial"), se=FALSE) +
xlab("Mean Chronic Conditions (Scaled)") +
ylab("Number of Visits") +
theme_classic()## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
ggplot(county_subset, aes(x = ELIX_SCALED, y = TOTAL_VISITS)) +
geom_point() +
stat_smooth(method = "lm", aes(color="linear"), se=FALSE) +
stat_smooth(method = "glm", aes(color="gamma"), method.args=list(family = Gamma(link = "log")), se=FALSE) +
stat_smooth(method = "glm.nb", aes(color="negative binomial"), se=FALSE) +
xlab("Avg Elix Score (Scaled)") +
ylab("Total Visits") +
theme_classic()## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
ggplot(county_subset, aes(x = LOS_SCALED, y = TOTAL_VISITS)) +
geom_point() +
stat_smooth(method = "lm", aes(color="linear"), se=FALSE) +
stat_smooth(method = "glm", aes(color="gamma"), method.args=list(family = Gamma(link = "log")), se=FALSE) +
stat_smooth(method = "glm.nb", aes(color="negative binomial"), se=FALSE) +
xlab("Length of Stay (Scaled)") +
ylab("Total Visits") +
theme_classic()## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace =
## control$trace > : iteration limit reached